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Abstract 

In  this  paper  we  evaluate  and  compare  two  algorithms  for  the  calculation  of  KAM  tori  in  Hamil¬ 
tonian  systems.  The  direct  fitting  of  a  torus  Fourier  series  to  a  numerically  integrated  trajectory 
is  the  first  method,  while  an  accelerated  finite  Fourier  transform  is  the  second  method.  The  finite 
Fourier  transform,  with  Hanning  window  functions,  is  by  far  superior  in  both  computational 
loading  and  numerical  accuracy.  Some  thoughts  on  applications  of  KAM  tori  are  offered. 

1.  Introduction 

The  KAM  theorem  has  laid  nearly  unnoticed  for  almost  sixty  years.  Named  for  Kolmogorov  [4],  Arnold 
[1],  and  Moser  [8],  this  theorem  predicts  that  most  lightly  perturbed  Hamiltonian  trajectories  should  lie  on 
the  surface  of  a  torus.  This  is  a  multiply-periodic  structure  that  should  be  represented  as  a  Fourier  series  in 
multiple  variables.  What  makes  this  different  from  a  perturbation  theory  is  that  the  Fourier  series  and  its 
basis  frequencies  need  not  be  approximate:  they  have  definite  values,  to  which  a  perturbation  theory  itself 
is  only  an  approximation. 

The  potential  applications  of  KAM  tori  are  numerous.  They  should  provide  the  general  solution  to  the 
problem  of  formation  flight  in  a  Hamiltonian  system.  They  can  be  used  as  another  way  to  compress  an 
ephemeris,  and  potentially  can  be  used  as  the  starting  point  for  perturbations  due  to  any  remaining  forces. 

To  achieve  this  promise,  we  must  have  an  efficient  way  to  produce  a  KAM  torus  passing  through  a  given 
set  of  initial  conditions,  and  assess  its  accuracy. 

2.  Two  Algorithms 

In  this  work  we  derive  and  compare  the  accuracy  and  performance  of  two  numerical  algorithms  for 
constructing  KAM  tori.  The  first  of  these  has  already  been  mentioned  in  the  literature,  and  consists  of 
numerical  fitting  a  KAM  torus  to  the  result  of  a  numerical  integration  of  an  orbit.  This  seems  to  be  the 
most  direct  method.  Second,  we  consider  extracting  the  Fourier  coefficients  of  the  torus  from  a  numerically 
integrated  Fourier  transform. 

Since  the  KAM  theorem  applies  directly  to  Hamiltonian  systems,  it  is  only  necessary  to  find  the  Fourier 
series  for  the  system  coordinates  g,;.  It  is  being  assumed  here  that  the  system  is  posed  in  a  set  of  “primitive” 
or  “native”  coordinates,  so  that  the  kinetic  energy  portion  of  the  Hamiltonian  is  quadratic  in  either  the 
generalized  velocities  qi,  or  equivalently  the  momenta  pi.  This  means  that  the  kinematic  half  of  Hamilton’s 
equations,  g.;  =  dTi/dpi  are  a  linear  set  of  equations 

q  =  A(q)p  +  b(q)  (1) 

in  the  momenta  p,  and  can  be  directly  solved  to  find 

p  =  A-1q  —  A-1b  (2) 

The  particular  forms  of  the  matrix  A{q)  and  the  vector  b (q)  depend  on  the  actual  Hamiltonian  system.  It 
is  to  be  noted,  however,  that  the  time  derivatives  of  the  coordinates  can  be  easily  extracted  from  the  same 
Fourier  series  that  produce  the  coordinates,  so  that  it  is  not  necessary  to  carry  separate  series  for  the  system 
momenta  pi. 

1The  views  expressed  are  those  of  the  author,  and  do  not  necessarily  reflect  those  of  the  United  States  Air  Force,  the 
Department  of  Defense,  or  the  U.S.  Government, 
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Also,  we  note  that  for  multiply  periodic  functions  f(t)  =  f(Q  =  oot)  that 
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is  a  norm  on  the  space  of  M  dimensional  multiply  periodic  functions.  It  satisfies  all  three  properties  required 
of  a  norm:  it  is  zero  if  and  only  if  /  =  0,  ||fc/(£)||  =  k  ||/(i)||  for  scalar  k,  and  ||  /(f)  +  g(t) ||  <  ||/(£)||  +  ||fl,(£)||- 
The  first  two  properties  are  trivial  to  show,  while  the  last  follows  easily  from  the  Fourier  series  form  of  the 
norm  using  the  equivalent  result  for  the  Euclidean  norm.  Note  that  invoking  the  ergodic  theorem,  Arnold 
and  Avez  [2] ,  the  space  integral  form  of  the  norm  (3)  can  be  evaluated  as  a  time  integral 
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for  Hamiltonian  systems.  We  will  find  this  result  useful  in  evaluating  how  successfully  we  have  approximated 
KAM  tori.  One  version  of  this  is  the  norm  of  the  error  in  the  equations  of  motion 
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where  the  torus  Fourier  series  supplies  the  coordinates  and  via  (2)  the  momenta,  and  a  time  derivative  of 
(2)  gives  the  predicted  p.  This  is  averaged  over  the  torus  via  (3). 

Finally,  the  notation  Q  =  ujt  is  to  indicate  that  these  variables  are  one  realization  of  the  Hamilton- Jacobi 
separation  coordinates. 


2.1  Direct  Fitting  of  Trajectories 


If  a  trajectory  lies  on  a  KAM  torus,  then  the  coordinates  can  be  represented  as  a  Fourier  series.  Since 
it  seems  the  most  obvious  thing  in  the  world  to  numerically  integrate  Hamilton’s  equations  to  produce  a 
trajectory,  fitting  the  result  of  such  a  numerical  integration  is  the  obvious  first  method.  It  has  been  done  in 
[7],  and  also  by  the  current  author,  [11]. 

Write  the  Fourier  series  representation  of  the  KAM  torus  as 


q(£)  =  C0  +  ^2  {Cj  cos(j  •  ut)  +  Sj  sin(j  •  ut)} 
j 


(6) 


The  multiple  summation  vector  jT  =  (ji,  j2,  ■■■Jn)  is,  of  course,  bounded  above,  and  the  KAM  torus  basis 
frequencies  are  most  conveniently  handled  as  a  vector  ujt  =  (oj\ .  oj2 ,  ...ojn)-  The  usual  rules  for  the  sine  and 
cosine  of  a  negative  argument  imply  that  the  first  nonzero  index  ji  should  be  positive  to  avoid  duplicate 
terms  in  the  series,  and  therefore  singular  least  squares  matrices.  After  the  “first”  positive  index,  subsequent 
indices  are  bounded  between  specified  truncation  limits  —jk,max  <  jk  <  jk,max ■  This  makes  for  a  somewhat 
cumbersome  index  summation  coding,  but  is  absolutely  essential  to  obtain  a  nonsingular  problem. 

In  the  direct  trajectory  fitting  method,  (6)  is  the  observable  relation  needed  for  least  squares  fitting.  Its 
linearization  is  then  easily  found,  the  matrix 


(7) 


where  the  unknowns  are  the  Fourier  coefficients  Cj,  Sj  as  well  as  the  system  basis  frequency  set  u>k-  The 
problem  is  nonlinear  with  the  inclusion  of  the  basis  frequency  set.  This  is  perhaps  the  most  general  case, 
since  each  coordinate  should  have  the  same  basis  frequencies,  and  that  forces  the  problem  into  the  nonlinear 
regime.  On  the  other  hand,  if  the  basis  frequencies  are  determined  separately  in  advance,  see  [5],  [11],  then 
each  coordinate  can  be  fit  separately,  and  the  algorithm  becomes  a  linear  least  squares  problem.  We  have 


not  adopted  this  option.  The  residuals  r.(  at  the  points  from  the  numerical  integration  are  then  calculated, 
and  corrections  to  the  state  are  given  by  the  usual  least  squares  formula 

5X  =  (TtT)  XTr  (8) 

Note  that  this  requires  the  inversion  of  ever  larger  matrices  as  the  order  of  the  Fourier  series  is  increased. 

2.2  Finite  Time  Fourier  Transforms 


The  usual  definition  of  the  Fourier  transform  of  a  function  of  time  f(t)  is 

/OO 

f(t)e~2ni,itdt 

-OO 


(9) 


In  any  numerical  application,  of  course,  the  use  of  infinite  time  intervals  is  out  of  the  question,  and  the 
Fourier  transform  of  a  multiply  periodic  function  becomes  an  inconvenient  sum  of  delta  functions,  which  are 
difficult  to  handle  numerically.  It  then  becomes  natural  to  investigate  the  finite  time  Fourier  transform 

Triy)  =  ^  f_Tme-^vitxP{t)dt 

=  JfjJ  f®6-2™***®*1*  ~JfJT  fWe~2nVitXrWdt  (10) 

While  the  Fourier  transform  is  defined  in  terms  of  the  cycle  frequency  v,  we  will  most  often  cite  results 
in  terms  of  the  angular  frequency  to  =  2/kv.  In  our  application  the  functions  f(t)  will  be  coordinates  of  a 
Hamiltonian  system  obtained  by  numerical  integration  of  a  trajectory.  The  symmetric  time  interval  form  of 
the  Fourier  transform  requires  two  separate  integration  intervals  when  starting  from  a  nominal  t  =  0.  While 
this  is  somewhat  of  a  complication,  it  does  make  best  use  of  the  finite  accuracy  available  through  a  numerical 
integration,  and  there  are  other  advantages  as  well. 

Following  Laskar  [6]  we  have  written  Xp{t)  f°r  a  cosine  window,  or  Hanning  window  function,  where 
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Laskar  has  studied  the  accelerated  convergence  of  the  frequency  determination  when  a  window  function  is 
used.  We  will  see  shortly  that  they  are  also  extremely  advantageous  in  decoupling  spectral  line  amplitudes. 
We  have  found  that  it  is  very  desirable  to  produce  the  numerical  integration  of  the  orbit  f(t)  as  a  first  step, 
and  then  to  import  it  into  the  Fourier  transform  code  as  a  very  large  vector  /(tj).  This  allows  efficient 
mechanization  of  (10)  by,  e.g.  Simpson’s  rule,  and  also  allows  multiple  analysis  options  without  repeating 
the  orbit  integration. 

Fig.  1  shows  the  Fourier  transform  of  a  single  spectral  line  using  the  first  four  Hanning  acceleration 
functions.  The  finite  Fourier  transform  oscillates  with  a  frequency  of 


lot  =  7t/T 


(12) 


so  this  plot  shows  only  the  immediate  vicinity  of  the  spectral  line  at  u>q.  There  are  several  things  to  note  from 
this  plot.  First,  the  peak  amplitude  is  unaltered  by  using  different  Hanning  functions.  Simple  calculation 
gives  the  Fourier  series  coefficients  for  a  single  line  as 

C  =  25L40,  S  =  2%A0  (13) 

where  Ao  is  the  peak  amplitude.  (The  constant  term  can  be  found  either  as  Co  =  :ftA(uJ  =  0),  or  from 
the  initial  conditions.)  Second,  while  it  is  well  known  that  higher  order  Hanning  functions  accelerate  the 
convergence  properties  of  finding  the  peak  frequency,  they  also  broaden  the  peak.  This  expands  the  realm  of 
convergence  to  the  actual  spectral  line  frequency  loq  without  confusion  by  sidelobes.  Finally,  recall  that  the 
entire  spectrum  will  be  an  infinite  summation  of  such  spectral  lines.  We  have  found  that  without  a  Hanning 
function,  the  amplitude  of  one  line  Ao  can  be  significantly  changed  by  the  far  tails  of  adjacent  lines.  This 


Figure  1: 


Fourier  transform  of  a  single  spectral  line  with  unit  amplitude,  using  different  Hanning  functions. 


makes  it  necessary  to  solve  a  large  system  of  linear  equations.  On  the  other  hand,  by  the  time  p  =  2,  which 
we  now  use  for  all  cases,  the  sidelobes  decay  very,  very  quickly.  This  means  that  we  can  directly  use  (13) 
for  all  but  the  very  closest  adjacent  lines.  (Parenthetically,  note  that  if  two  lines  are  very  close,  so  that 
j  ■  ss  k  ■  <*>,  then 

(j  -  k)  •  «  «  0  (14) 

which  is  just  another  form  of  the  resonance  condition.  So,  if  a  KAM  torus  exists  (no  resonance),  the  above 
will  not  be  a  problem.)  This  is  a  huge  advantage  of  the  finite  Fourier  transform  method,  since  the  decoupling 
of  spectral  lines  makes  it  possible  to  solve  very  large  problems  one  spectral  line  at  a  time. 

Finally  another  advantage  of  the  finite  Fourier  transform  method  is  that  it  does  not  require  an  approxi¬ 
mation  to  the  actual  KAM  Fourier  series,  as  does  the  nonlinear  trajectory  fitting  method.  Rather,  only  an 
approximate  knowledge  of  the  system  frequencies  is  needed. 

3.  Numerical  Results 

In  this  effort  we  have  studied  a  series  of  tori  in  the  restricted  problem  of  three  bodies,  see  e.g.  [10].  The 
Hamiltonian  function  is  given  by 

n  =  \  (pI  +  pI)  +  vpx  -  xpy  -  ^  (15) 

in  the  usual  rotating  frame  of  reference.  Here  p  is  the  mass  of  the  smaller  primary,  and  the  radii  are  given 

by 

r\  =  (x  -  pf  +  y2,  r\  =  (a;  +  1  -  pf  +  y2  (16) 

One  advantage  of  studying  this  system  is  that  we  can  use  the  surface  of  section  technique  to  visualize  the  tori, 
[3],  We  have  studied  a  group  of  tori  for  a  value  of  H  =  —1.6  with  p  =  0.01214,  nearly  the  earth  /  moon  ratio, 
and  the  usual  dimensionless  units  being  used.  The  periodic  orbit  at  the  core  of  the  tori  has  initial  conditions 


Figure  2:  A  three  dimensional  view  of  a  KAM  torus. 


x  =  0.55954260514673,  py  =  1.4186361935797  with  y  =  px  =  0,  and  a  period  of  r  =  5.597166681940  time 
units  (TU).  This  gives  the  single  torus  frequency  for  the  periodic  orbit  as  u \  =  1.1225653378258  radians  per 
time  unit.  While  the  periodic  orbit  itself  has  only  one  frequency,  nearby  tori  should  approximate  the  Floquet 
solution,  so  the  second  frequency  should  be  nearly  the  Poincare  exponent  up  =  0.1592640457  radians  /  TU 
for  this  orbit.  Of  course,  the  Fourier  series  for  the  periodic  orbit  can  be  obtained  by  standard  numerical 
methods.  Finally,  we  note  that  this  group  of  tori  is  part  of  the  2:1  resonance  structure  for  this  mass  value. 

Most  of  the  results  that  follow  are  for  the  torus  that  has  initial  conditions  x  =  0.600,  y  =  0,  px  =  0, 
py  =  1.3322289632022.  A  view  of  the  final  constructed  torus  is  shown  in  Fig.  2.  Blue  and  green  lines  are 
lines  of  constant  Q\  and  Q2,  while  the  trajectory  is  shown  in  red.  While  this  can  appear  complicated  when 
projected  back  into  the  physical  phase  space,  the  transformation  actually  makes  the  dynamics  quite  simple. 
All  the  Qi  coordinates  are  linear  in  time.  Fig.  3  shows  the  behavior  in  the  torus  coordinate  space.  In  fact, 
this  is  the  “surface  texture”  used  in  the  three  dimensional  model  of  Fig.  2 

One  point  that  needs  to  be  brought  out  here  is  that  all  methods  depend  on  knowing  the  torus  frequencies 
with  relatively  good  accuracy  when  beginning  the  fit,  by  any  of  the  methods  discussed.  These  use  relatively 
large  final  times  tmax  for  a  numerical  integration  in  order  to  sample  the  torus  are  required  when  attempting 
to  sample  uniformly  across  the  Q  space.  While  no  numerical  integration  is  required  fitting  the  equations  of 
motion.  In  order  for  the  linearization  we  have  used  to  be  valid,  this  means  that  the  error  in  the  separation 
coordinates  5Qi  ss  5uitmax  «  1.  Since  the  final  time  may  be  relatively  large,  this  implies  that  the  frequency 
error  must  be  correspondingly  small  when  even  beginning  the  process.  The  author  suspects  that  fits  that 
do  not  converge  convincingly  to  small  errors  may  have  “hung”  with  an  error  in  some  phases  5Qi  that  is  not 
small. 

The  question  of  how  far  to  extend  the  Fourier  solution  naturally  arises.  Fig.  4  shows  the  result  of  fitting 
a  restricted  problem  torus  using  the  least  squares  process,  and  by  the  finite  Fourier  transform.  The  orbit 
fit  solution  uniformly  improves  until  a  limit  is  reached,  beyond  which  improvement  comes  more  slowly,  if  at 
all.  This  may  be  associated  with  the  accuracy  limits  involved  in  the  solution  of  a  very  large  linear  system  of 
equations  in  the  case  of  the  least  squares  orbit  fit.  The  finite  Fourier  transform  results  continue  to  improve 
for  some  considerable  time  after  the  orbit  fitting  process  has  reached  an  accuracy  floor.  Since  the  Fourier 
transform  computational  burden  grows  only  linearly  with  the  number  of  spectral  lines  included,  it  has  a 
clear  advantage  here  both  in  terms  of  accuracy  and  in  terms  of  computational  loading. 

The  residuals  in  the  orbit  points  assumes  an  interesting  form  once  the  match  between  the  torus  and  the 
numerical  integration  is  close  enough.  When  the  two  are  still  significantly  apart,  definite  patterns  appear  in 
the  residuals,  and  the  power  spectral  density  plot  of  the  residuals  will  often  reveal  where  significant  structure 
is  yet  to  be  accounted  for.  But  when  close,  the  coordinate  residuals  often  assume  a  form  as  shown  in  Fig 
5.  This  linear  trend  can  only  be  accounted  for  by  an  error  in  one  or  more  fundamental  torus  frequencies.  A 
missing  harmonic  would,  of  course,  produce  a  periodic  signature  in  the  residuals.  But  an  error  in  a  frequency, 
if  small  enough,  shows  up  as 


Sk  sin((wo  +  6u)t) 


Sk  sin(o;ot)  cos (Sut)  +  Sk  cos(wot)  sin(i5u;f) 


Figure  3:  Motion  on  the  torus  surface  is  simple  straight  line  motion. 


Finite  Fourier  _  Orbit  Fit 


Figure  4:  Equation  of  motion  error  norm  for  direct  fitting  of  a  torus  from  a  numerical  integration.  Red  is 
the  Fourier  transform,  while  Green  shows  the  orbit  fit  results. 


~  Sk  sin(wot)  +  Sk  cos(u>ot)5ujt  (17) 

This  produces  a  secular  term  in  the  residuals.  The  error  shown  in  Fig  5,  which  has  a  growth  of  a  few  parts 
in  106  over  a  time  interval  of  several  thousand  time  units  implies  an  error  in  a  frequency  of  a  few  parts  in 
109.  This  may  be  in  the  torus  decomposition,  or  an  error  in  the  numerical  integration.  But  it  is  a  very  small 
effect.  And  in  fact,  it  may  arise  from  the  numerical  integration,  since  the  discrete  jumps  in  the  slope  cannot 
come  from  the  Fourier  series. 


4.  Action  Momenta 

The  Fourier  series  angles  Qi  =  flit  are  the  N  separation  coordinates  expected  from  the  Hamilton- Jacobi 
theorem.  Their  conjugate  momenta  are  the  action  variables  [9] 

T>i=  ^  <j>  p  •  dq  (18) 

where  the  integral  is  carried  out  around  the  N  irreducible  contours  Tj  around  the  torus.  This  is  easily  done 
with  the  torus  Fourier  series  in  hand,  since  we  integrate  around  one  new  coordinate  Qi,  holding  all  the  others 
constant.  It  is  also  known  that  these  integrals  are  invariant  to  the  particular  realization  of  these  contours. 
Starting  from  the  fitted  Fourier  series  for  the  coordinate  vector  q,  the  system  momenta  can  be  found  through 
the  kinematic  half  of  Hamilton’s  equations,  q  =  dTi/dp.  This  often  can  be  rearranged  into  the  form 

p  =  T2q  +  Tiq  (19) 

where  the  matrix  T2  comes  from  the  quadratic  part  of  the  kinetic  energy,  and  the  matrix  Tj  most  often 
appears  when  a  rotating  reference  frame  is  employed.  Of  course,  this  is  not  the  most  general  form  of  the 
system  momenta,  but  will  be  sufficient  for  this  paper.  Then  the  coordinate  rates  are  obtained  by  a  time 
derivative 

q  =  Ci '  w)  {~CJ  sin0 ' 9)  +  Sj  cos(j  •  0)}  (20) 

j 

where  we  have  extended  the  scalar  definition  (6)  to  a  vector  of  coordinates  q.  Also,  then,  we  can  write  the 
differential  dq  for  contour  T.;  as 
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Figure  5:  Coordinate  residuals  when  the  torus  fitting  error  is  dominated  by  a  frequency  error  and  frequency 
jumps. 


(21) 


=  •?’*  {  _Cj  Sill(j  '  &)  +  SJ  C0S(j  '  3)  }  ddi 

j 

So,  having  fit  the  KAM  coordinates  q  and  having  obtained  the  basis  frequencies  fij,  it  is  then  possible  to 
numerically  calculate  the  conjugate  action  momenta  V, . 

4.  Conclusions  and  Prospectus 

In  this  paper  it  has  been  shown  that  the  finite  Fourier  transform  algorithm  has  very  great  advantages 
in  both  accuracy  and  computational  loading  for  reducing  numerically  integrated  orbits  to  KAM  Tori.  This 
process  has  been  carried  out  to  better  than  single  precision  accuracy,  and  the  resulting  Fourier  series  allow 
the  reduction  of  the  dynamics  to  local  Hamilton- Jacobi  form:  the  coordinates  increment  linearly  with  time, 
while  the  momenta  are  constant. 

Research  currently  is  proceeding  in  applying  KAM  torus  theory  to  earth  satellites,  including  the  prob¬ 
lems  of  formation  flight,  stationkeeping  satellite  constellations,  compressing  ephemeris  data  for  navigational 
satellites,  and  the  important  problem  of  determining  a  torus  from  observational  tracking  data.  This  latter 
problem  involves  the  solution  to  the  problem  of  motion  in  the  vicinity  of  a  KAM  torus,  which  the  author 
hopes  to  report  on  shortly. 
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